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Abstract 

Thermodynamic perturbation theory is applied to the model of globular proteins studied by 
ten Wolde and Frenkel (Science 277, pg. 1976) using computer simulation. It is found that the 
reported phase diagrams are accurately reproduced. The calculations show how the phase diagram 
can be tuned as a function of the lengthscale of the potential. 
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I. INTRODUCTION 



One of the most important problems in biophysics is the characterization of the struc- 
ture of proteins. The experimental determination of protein structure by means of x-ray 
diffraction requires that the proteins be prepared as good quality crystals which turns out 
to be difficult to achieve. Given the fact that recent years have seen an explosion in the 
number of proteins which have been isolated, the need is therefore greater than ever for 
efficient methods to produce such crystals. Without finely tuned experimental conditions, 
often discovered through laborious trial and error, crystallization may not occur on labo- 
ratory time-scales or amorphous, rather than crystalline, structures may form. The recent 
observation by ten Wolde and Frenkell] of enhanced nucleation of a model protein in the 
vicinity of a metastable critical point is thus of great interest and could lead to more efficient 
means of crystallization if such conditions can be easily identified for a given protein. 

Wilson noted that favorable conditions for crystallization are correlated with the behavior 
of the osmotic second virial coefficient [2f and, hence, depend sensitively on temperature. If 
the second virial coefficient is too large, crystallization occurs slowly and if it is too small, 
amorphous solids form. By comparing the experimentally determined precipitation bound- 
aries for several different globular proteins as a function of interaction range, controlled 
by means of the background ionic strength, Rosenbaum et al have shown that the phase 
diagrams of a large class of globular proteins can be mapped onto those of simple fluids 
interacting via central force potentials consisting of hard cores and short-ranged attractive 
tails . They also discuss the important fact that the range of interaction can be tuned 
by varying the composition of the solvents used. The attraction must, in general, be very 
short ranged if this model is to apply since a fluid-fluid phase transition is not typically 
observe experimental and it is know,, that this transition is on.y suppressed in simple 
fluids when the attractions are very short ranged 5] . These studies therefore support the 
conclusion that the study of simple fluids interacting via potentials with short-ranged attrac- 
tive tails can give insight into nucleation of the crystalline phase of a large class of globular 
proteins, ten Wolde and Frenkel have studied nucleation of a particular model globular 
protein consisting of a hard-core and a modified Lennard- Jones tail by direct free energy 
measurements obtained from computer simulations[l]. They found that the nucleation rate 
of a stable FCC solid phase could be significantly enhanced in the vicinity of a metastable 
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critical point. The enhancement is due to the possibility that a density fluctuation in the 
vapor phase is able to first nucleate a metastable droplet of denser fluid which, in turn, 
forms a crystal nucleus. The fact that intermediate metastable states can accelerate barrier 
crossing has been confirmed using kinetic models :&] and the physics of the proposed non- 
classical nucleation model has also been confirmed by theoretical studies based on density 
functional models[3],j^. This observation opens up the possibility of efficiently producing 
good quality protein crystals provided that it is understood how to tune the interactions 
governing a given protein so that its phase diagram possesses such a metastable state under 
experimentally accessible conditions. A prerequisite for achieving this is to go beyond the 
heavily parameterized studies conducted so far and to be able to accurately predict phase 
diagrams given knowledge of the range of the protein interactions. In this paper, we describe 
the application of thermodynamic perturbation theory to calculate the phase diagram based 
solely on the interaction model. In so far as the range of interaction is important, and not 
the detailed functional forms, this approach, if successful, gives a direct connection between 
the phase diagram and the range of interaction without the need for further, phenomenolog- 
ical parameterizations. We show that the theory can be used to successfully reproduce the 
phase diagrams of ten Wolde and Frenkel based only on the interaction potential and assess 
the effect of the range of the interatomic potential on the structure of the phase diagram. 

In the next Section, the formalism used in our calculations is outlined. This involves the 
standard Weeks- Chandler- Andersen perturbation theory with modifications to improve its 
accuracy at high densities. The third Section discusses the application of the perturbation 
theory to the ten Wolde-Frenkel interaction model. Whether or not perturbation theory 
is applicable to this type of system is not immediately evident: the hard-core square well 
potential has long served as a test case for developments in perturbation theory 9]. So we 
show how the size of the various contributions to the total free energy varies with temperature 
and that second order contributions to the free energy are of negligible importance. In the 
fourth Section, the calculated phase diagram for the hard core plus modified Lennard- Jones 
tail is shown to be in good agreement with the reported Monte Carlo (MC) results. Since 
the perturbation theory is also well known jlo| to give a good description of long-ranged 
potentials such as the standard Lennard- Jones, we expect that it can be used with some 
confidence to explore the effect of the length scale of the potential on the phase behavior 
of the systems. To illustrate, we present the phase diagram as a function of the range 
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of the modified Lennard- Jones tail and show that the appearance of the metastable state 
requires only a minor modification of the range of the potential. The final Section contains 
our conclusions where we discuss the prospect for using the perturbation theory free energy 
function as the basis for density functional studies of the nucleation process and for studies 
of the effect of fluctuations on the nucleation rate. 



II. THERMODYNAMIC PERTURBATION THEORY 

Thermodynamic perturbation theory allows one to express the Helmholtz free energy F 
of a system in terms of a perturbative expansion about some reference system. There are 
a number of different approaches to constructing the pert urbative expansions such as the 
well known Weeks- Chandler- Andersen ( WCA) jll | . jl^] . . [l^ theory and the more recent 
Mansoori-Canfield/Rasaiah-Stell theory [lsj]. The latter appears to be more accurate for 
systems with soft repulsions at small separations while the former works better for systems 
with stronger repulsions. Here, we will be interested in a hard-core potential with a modified 
Lennard- Jones tail, so we use the WCA theory as modified by Ree and coworkers . jlo|] . ~v\ 
as discussed below. The first step is to divide the potential into a (mostly) repulsive short- 
ranged part and a longer ranged (mostly) attractive tail according to the prescription 

v(r) = v (r)+w(r) (1) 
v (r) — v (r ) — v' (r ) (r — r ) , r < r 
0, r > r 

v (r ) + v' (r ) (r - r ) , r < r 
v (r) , r > r 



v (r) 



w (r) 



The short ranged part is generally repulsive and can therefore be well approximated by a 
hard-sphere reference system. The long-ranged tail describes the attractive forces and must 
also be accounted for so that distinct liquid and gas phases exist (i.e. so that the phase 
diagram exhibits a Van der Waals loop). There are a number of versions of the WCA- 
type perturbation theory depending on the choice of the separation point r . Barker and 
Henderson chose the separation point tq to be the point at which the potential goes 
to zero, v (tq) = 0, (they also did not include the linear term in the expressions above). 
Subsequently, WCA achieved a better description of the Lennard- Jones phase diagram by 
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taking the separation point to be at the minimum of the potential, v' (tq) = . Reefl^ first 
suggested that the free energy be minimized with respect to ro, and introduced the linear 
terms in eq. (0) , in which case the first-order perturbation theory is equivalent to a variational 
theory based on the Gibbs-Bugolyubov inequalities wj]. Later, Ree and coworkers showed 
that essentially the same results could be achieved with the prescription 

r = min (r min , r nn ) (2) 

where r min is the minimum of the potential, v' (r^n) = and r nn = 2sp -1 / 3 , where p is 
the density, is the FCC nearest-neighbor distance|ll|. For low densities, this amounts to 
the original WCA prescription whereas for higher densities, the separation point decreases 
with increasing density. In this case, the linear term in the definition of Vq (r) ensures 
the continuity of the first derivative of the potential. Calculations for the Lennard- Jones 
potential, as well as inverse power potentials, show that this modification of the original 
WCA theory gives improved results at high density. Finally, eq.(j2} was modified to switch 
smoothly from r min to r nn as the density increases so as to avoid discontinuities in the free 
energy as a function of density and thus singularities in the pressure [3]. Below, we will 
refer to this final form of the Weeks-Chandler-Andersen-Ree theory as the WCAR theory. 

A. Contribution of the long-ranged part or the potential 

The contribution of the long-ranged part of the potential to the free energy is handled 
perturbatively in the so-called high-temperature expansion 

jfPF - I/JF. = i/J (W)„ + ±f \{W\ - (W)l\ + ... (3) 

where F Q is the free energy of a system of N particles subject only to the short-ranged 
potential Vq (r) at inverse temperature (3 = and where the total attractive energy is 

w= E w (ra) ■ ( 4 ) 

l<i<j<N 

The brackets (...) indicate an equilibrium average over a system interacting with the po- 
tential Vq. The first term on the right is easily calculated since it only involves the pair 
distribution function of the reference system 

(W) = \p P J dt g (r) w (r) (5) 



where go (r) is the pair distribution function of the reference system. The second term 
requires knowledge of three- and four-body correlations for which good approximations are 
not available. Its value is typically estimated using Barker and Henderson's "macroscopic 
compressibility" approximation [ 
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" Wo] * -\P 2 P i-^p) J g (r) w 2 (r) (6) 
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where Pq is the pressure of the reference system at temperature ksT = 1/(3 and density p 



B. Contribution of the short-ranged part of the potential 

The description of the reference system is again accomplished by perturbation theory. 
Since the potential vq (r) is not very different from a hard core potential, this perturbation 
theory does not involve the high temperature expansion but, rather, involves a functional 
expansion in the quantity exp (— (5vo (r)) — exp (— /3vh s (r; d)) where Vh s (r;d) is the hard 
sphere potential for a hard-sphere diameter d. The result is 

^(3Fo - ^(3F hs (pd 3 ) = j dt y hs (r) (exp {-f3vo (r)) - exp {-(3v hs (r; d))) + ... (7) 

where yh s (r, pd 3 ) is the hard-sphere cavity function, related to the pair distribution function 
through ghs( r ) = exp (— (3vh s ( r ; d)) y^ s (r). Several methods of choosing the hard-sphere 
diameter of the reference system are common. The WCA prescription is to force the first 
order term to vanish 

= J d~t y hs (r) (exp {-(3v (r)) - exp (-/3v hs (r; d))) . (8) 

and a simple expansion about r = d gives the cruder Barker and Henderson 
approximation which gives 

J dr (exp (-Pv (r)) - 1) + J dr (1 - exp {-(3v hs (r; d))) ~ 0. (9) 

In either case, one can then consistently approximate the pair distribution function of the 
reference state as either 

9o (r) ^ g hs (r) (10) 

or 

g (r) ~ exp(-(3v (r))y hs (r) (11) 



where the difference between using one expression or the other is of the same size as ne- 
glected terms in the perturbation theory. Here, we follow Ree et al[li3] in using the WCA 
prescription for the hard-sphere diameter, eq.(jSJ) and the first approximation, eq. (jlU)) . for 
the pair distribution function. Then, the complete expression for the free energy becomes 

^(3F = ^3F hs (pd 3 ) + J ' d~f y hs (r) (exp (-8v (r)) - exp (-(3v hs (r; d))) (12) 
+\@P J 9hs (r)w(r) 
~\fiP \ ^7T I / dy¥ 9hs (r) w 2 (r) . 
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The pressure, P, and chemical potential p are calculated from the free energy using the 
standard thermodynamic relations 

p opN 

1 BP 

Pp = T7^+— • 
N p 

C. Description of the reference liquid 

The calculation of liquid phase free energies require as input the properties of the hard 
sphere liquid. These are known to a high degree of accuracy and introduce no significant 
uncertainty, nor any new parameters, into the calculations. 

The properties of low density hard-sphere liquids are well described by the Percus-Yevick 
(PY) approximation but this is not adequate for the dense liquids to be considered here. So 
for the hard sphere cavity function, we have used the model of Henderson and GrundkejsjJ 
which modifies the PY description so as to more accurately describe dense liquids. The 
corresponding pair distribution function is then that of Verlet and Weiss and the equation 
of state, as obtained from it by both the compressibility equation and the pressure equation, 
is the Carnahan-Starling equation of state[l4|. The free energy as a function of density 
follows immediately and is given by 

1 or, f.J*\_^f^*\ 1 , _ 4-377 



-p Fhs ^)=ln( P A*)-l + V ^—^ (14) 

where rj = pd 3 . The second term of eq.([12|) is easily evaluated numerically as its kernel is 
sharply peaked about r = d . The most troublesome part of the calculation is the evaluation 
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of the contributions of the long-ranged part of the potential, w (r) . One method is to divide 
the necessary integrals along the lines 

J d?? g hs (r) w (r) = J d~r* w (r) + J d~r* (g hs (r) - 1) w (r) (15) 

where the first piece can be calculated analytically and the second involves the structure 
function g^ s (r)—l which is relatively short ranged allowing a numerical evaluation. However, 
at high densities this can still be difficult to evaluate as the hard-sphere structure extends for 
considerable distances. In the Appendix, we discuss a more efficient method of evaluation 
based on Laplace transform techniques. We have used both methods and obtained consistent 
results: in general, the second is much easier to implement and numerically more stable. 



D. Description of the reference solid 

To calculate the properties of the solid phase, the same expressions are used except that 
the reference free energy is now that of the hard-sphere solid and the pair distribution 
function is the spherical average of the hard-sphere pair distribution function. Both of these 
quantities can be obtained by means of classical density functional theory, but here we 
choose the simpler, and older, approach which makes use of analytic fits to the results of 
computer simulations together with the known high-density limit of the equation of state. 
This limits the present calculations to the investigation of the FCC solid phase as this is 
the only one for which extensive simulations have been performed. We stress that these fits 
are very good and that they introduce no new parameters into the calculations of the phase 
diagrams. 

In the calculations presented below, we have used the equation of state proposed by 
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— = 3—^— + 2.557696 + 0.12530776 + 0.17623936 2 - 1.0533086 3 (16) 
P Vc-V 

+2.8186216 4 - 2.9219346 5 + 1.1184136 6 

where r) c = fx/2 is the value of the packing fraction at close packing. Notice that the first 
term is the high density limit of the Lennard- Jones-Devonshire cell theory which is expected 
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to be exact near close packing (see, e.g., the discussion in 23]). The free energy is then 
calculated by integrating from the desired density to the close-packing density giving 



(3F = (3F 



P 



LJDA 



dl] 

V ' 



(17) 



For the spherically-averaged pair distribution function for the FCC solid, we use the 



analytic fits of Kincaid and Weiss 



24] 



9kw ( r ) = {A/x) exp — w\ (x — x±) — w\ (x — xif 



+ 



w 



E 



11; 



i=2 XiX 



oxp (— w 2 (x — X\f^j 



(18) 

Here x = r/d, the parameter A is fixed by requiring that the pressure equation reproduce 
the Hall equation of state 

(19) 



— J = 1 + 4r}g KW {!) 



P 



Hall 



and the parameters Wi, w 2 and w are given as functions of density by analytic fits to the 



MC data[24j. No such fit is given for the parameter x± so its value must be determined by 
interpolating from the values extracted from the MC data as given in [24 1 . The quantities 
rij and Xi are the number of neighbors and the position of the z-th lattice shell respectively. 
Note that Ree et al suggest using the earlier parameterization of Weisj^] at lower densities, 
where it is slightly more accurate, and the Kincaid- Weis version at higher densities. We have 
not done this because it leads to discontinuities in the free energy as a function of density 
at the point the switch is made. Since these are just empirical fits, we do not believe there 
is a significant loss of accuracy. 



III. APPLICATION TO A MODEL PROTEIN INTERMOLECULAR POTEN- 
TIAL 

A. The potential 

The only input needed for the perturbative calculation outlined in Section II is the in- 
termolecular potential: there are no phenomenological parameters to specify. The goal of 
this work is to show how to construct a realistic free energy functional with which to study 
nucleation of protein crystallization using the model potential of ten Wolde and Frenkel[l|. 
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This interaction model consists of a hard-sphere pair potential with an attractive tail 

oo , r < a 



v (r) 



4e 



"ff - "7 \ \? 



(20) 

, r > a 



((f) 2 - 1 ) ((J) 2 - 1 ) 

The tail is actually a modified Lennard- Jones potential and the two are related by 

G(r-a)v(r)=e(r- a) v LJ ^a^-j - lj . (21) 

As such, the potential decays as a power law and is not short-ranged in the usual sense. 
Nevertheless, as a becomes larger, the range of the potential decreases: for example, the 
minimum of the potential is 



1 + ( t ) (22) 



o V \aj 
v(r min ) = -e. 

which approaches the hard core for large a . Furthermore, for a fixed position r > r min , the 
value of the potential decreases with increasing a relative to its minimum. For example, 

v(2a)/v(r min )= — 2 (23) 

so that as a increases, the interactions of particles separated by much more than r min con- 
tribute less and less to the total energy compared to the contribution of particles that are 
close together. Figure 1 shows the evolution of the shape of the potential as a increases. 
The range of the potential varies from about 2.5 hard sphere diameters for a — 1 to less 
than 1.25 diameters for a = 50. Also shown in the figure is the separation of the potential 
into long- and short-ranged pieces for the case a = 50 where it is clear that even for this very 
short-ranged potential, the long-ranged function W(r) varies relatively slowly compared to 
the short-ranged repulsive potential Vo (r). 

B. Comparison of various approximations 

Figure 2 shows the contributions of the various terms contributing to the free energy 
at two temperatures. In both cases, the second order term is seen to be negligible. This 
is because at low density, the free energy is dominated by the ideal-gas contribution, all 
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FIG. 1: The ten Wolde-Frekel potential as a function of a. The inset shows the division of the 
potential into long-ranged and short-ranged parts for a = 50. 

other contributions going to zero with the density, whereas at moderate to high density, the 
compressibility controlling the size of the contribution of the second order term, see eq.([12p. 
diminishes quickly from its zero density limit of 1.0 to something on the order of 0.1 at 
moderate densities and is of order 0.01 at high densities. We conclude that the second order 
contributions, at least calculated within the macroscopic compressibility approximation, 
eq.(jnj), can be neglected. 

In the case of the lower temperature, fcgT '/e = 0.35, the first order contributions quickly 
grow with density until at high densities, they are larger than the zeroth order contributions 
thus suggesting that the perturbation theory will not prove very accurate. At the higher 
temperature, /c#T ' je = 1.5, the first order contributions are much better controlled and we 
expect the perturbation theory to be relatively accurate. 

We have also tested the various approaches to the selection of the separation point of the 
potential - the WCA prescription, eq.(J2J), the WCAR prescription and minimization of the 
free energy with respect to r . As expected, the only significant differences occur at high 
density, where variations of the free energy of 10% occur, but we find virtually no effect on 
the phase diagram. 
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FIG. 2: The various terms contributing to the total free energy as a function of density for two 
different temperatures. At the lower temperature, the first order contribution dominates the hard- 
sphere contribution whereas at higher temperatures, the zeroth order terms dominate. 

IV. PHASE DIAGRAMS 

Figure 3 shows the phase diagram as calculated from the WCAR theory for the potential 
and parameters used by ten Wolde and Frenkel[l[ and its comparison to the results of 
Monte Carlo simulations by these authors for a = 50. The lines, from our calculations, 
and the symbols, from the simulations, divide the density-temperature phase diagram into 
three parts: the liquid region (low density and high temperature), the fluid-solid coexistence 
region and the solid region (at high density). In the calculations, the lines are determined by 
finding, for a given temperature, the liquid and solid densities that give equal pressures and 
chemical potentials for the two phases as determined using eq. f|13[) based on the liquid and 
solid free energy calculations (which differ only in the equation of state and pair distribution 
function of the reference states). The fluid-fluid coexistence is determined similarly except 
that the free energy for both phases is calculated using the same reference state (the hard- 
sphere fluid) with the resulting free energy exhibiting a Van der Waals loop. 

The calculations and simulations are in good qualitative agreement with a fluid-fluid 
critical point that is suppressed by the fluid-solid phase boundaries. The values of the 



12 




FIG. 3: Comparison of the predicted phase diagram, lines, to the Monte Carlo results, symbols, of 
ref.Q for a = 50. Some error bars are superposed on the symbols. 

coexisting densities are in good agreement at low temperatures, where the liquid density 
is very low and at high temperatures. That these limits agree is as expected from our 
discussion of the relative sizes of the various contributions to the free energy. It is perhaps 
surprising that the agreement is so good even for temperatures as low as ksT/e ~ 1. The 
intermediate temperature values, where the attractive tail and finite density effects are 
important, are the most poorly described. The same is true of the fluid-fluid coexistence 
curve. The critical point is estimated to occur at about ksT/e ~ 0.48 and pa 3 ~ 0.4 whereas 
the simulation results are ksT/e ~ 0.4 and pa 3 ~ 0.3. We have tested these results by using 
different choices for the pair distribution function of the reference state (see eqs. (ll0|) -(fTT j) ). 
and different choices for the division of the potential (such as minimizing the free energy 
with respect to the break point) but none of these alternatives produces any significant 
change. 

An interesting feature of short-ranged interactions is that under some circumstances, they 
give rise to solid-solid transitions where the lattice structure remains the same but solids of 
different densities can coexist (i.e. a van der Waals loop occurs in the solid free energy) 26^ . 
We have searched for, but find no evidence of, such a transition with the present potential. 

To give some idea of the typical energy barrier between the coexisting phases, we show in 
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FIG. 4: Calculated free energies as a function of density for the liquid and solid phases at ksT/e = 
0.4 and for a = 50. The points mark the location of the coexisting phases. 

Fig. 0]the calculated isothermal free energies as a function of density between the coexisting 
fluid and solid phases at fceT '/e = 0.4 for the short-ranged (a = 50) potential. The fluid has 
a density of 0.008 and Helmholtz free energy of -5.82 in reduced units. The maximum free 
energy is -2.57 and the solid free energy is -5.02 at a density of 0.88. 

Figure El shows the phase diagrams calculated from the WCAR theory as a function of 
the range of the potential (i.e., different values of a). For a = 1, for which the minimum of 
the potential well is r min — 1.5 and corresponding to a tail that closely resembles a standard 
Lennard- Jones interaction, the phase diagram has the classical form exhibiting three stable 
phases, a critical point and a triple point. As a increases, and the range of the potential 
decreases, the critical point moves towards the triple point. Even for a = 5 and r m i n = 1.31, 
the critical point lies very near the triple point and the two become nearly identical for 
a = 10 and r m i n = 1.26. Our conclusion is that for this model, the suppression of the triple 
point occurs when the range of the potential, as characterized by its minimum, falls to about 
a quarter of the hard-core diameter. 
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FIG. 5: Calculated phase diagrams as a function of a showing that the critical point is suppressed 
for a > 10. 

V. CONCLUSIONS 

Our aim here has been to provide a fundamental model of protein crystallization without 
the need for parameterizations other than the interaction potential. Since the potential 
for globular proteins can be tuned, by varying e.g. the background ionic strength of the 
solutions, this provides a rather direct connection between theoretical indications of favorable 
conditions for nucleation and experimentally accessible control parameters. 

We have shown that thermodynamic perturbation theory gives a good, semi-quantitative 
estimate of the phase diagram of a model interaction for globular proteins. The accuracy of 
the perturbation theory is expected to improve as the range of the potential increases so, e.g., 
the prediction of the value of a at which the critical point becomes suppressed is expected to 
be reasonably accurate. Unlike the results of a recent study of colloids interacting via short- 



ranged potentials 27|, we do not find that the second order terms in the high-temperature 
expansion play an important role in the structure of the phase diagram. 

This free energy calculation, which only uses the interaction model as input, should be 
contrasted with other more phenomenological approaches. In phase field models, the free 
energy is taken to be a function of one or more order parameters. The actual form of the 
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free energy is typically of the Landau form which is to say, a square-gradient term plus an 
algebraic function of more than second order in the order parameter. The coefficients of 
these terms must be fitted to experimental data and the adequacy of the assumed function 
is difficult to assess. Similarly, the recent density functional models of TalanquerfT] and 
Shiryayevpl] depend on an ad hoc free energy functional, based on the van der Waals free 
energy model for the fluid, with several phenomenological parameters. 

We believe that our work can serve as the basis for further theoretical study of the 
nucleation of globular proteins using density functional theory. While the present description 
of the two phases requires as input separate equations of state and pair distribution functions 
for the reference hard sphere fluid and solid phases, standard methods exist for interpolating 
between these so as to provide a single, unified free energy functional suitable to the study 
of free energy barriers (see,e.g. ref.|2ja|). Such a unified model can be used to study static 
properties, such as the structure of the critical nucleas, using density functional theory as 
well as the effect of fluctuations on the transition rates by the addition of noise obeying the 
fluctuation-dissipation theorem. 

Finally, it would be desirable to confront the approach developed here to experiments 
aiming to determine the interaction potential and the phase diagram of concrete globular 
proteins of interest such as lysozyme and catalase. In recent years, considerable effort 
was devoted to protein crystallization under microgravity conditions on the grounds that 
some undesirable effects such as densit y g radients and advection present in earth-bound 
experiments can be virtually suppressed [23]. In parallel, earth-bound experiments are being 
carried out to determine conditions and parameters to be used in a microgravity experiment. 
In either case, the role of the metastable critical point has so far not been addressed in a 
detailed manner. We believe that the availability of a theory as parameter-free as possible 
like the one developed in the present work could provide the frame for undertaking such a 
study on a rational basis. 
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APPENDIX A: EVALUATION OF LONG-RANGED CONTRIBUTION TO THE 
FREE ENERGY 

We begin by writing the first order contribution of the long-ranged potential as 

/ d~t g hs (r) w(r) = j d~t g hs (r) v (r) - J d~t g hs (r) v (r) (Al) 

so that the second term involves the very short ranged function v o (r) and is easily performed 
numerically. Our focus is therefore on the evaluation of the first term on the right. If we 
write the potential as the sum of a hard-core and a continuous tail 

v (r) = v hs (r) + 6 (r - a) v tail (r) (A2) 

and the effective hard-sphere diameter d > a, as it clearly will always be, then 

/ d~J* g hs (r) v (r) = d~r> (r - d) y hs (r) v (r) (A3) 

= J d~r> g hs (r) v taU (r) 

so that we can ignore the discontinuity of the hard-core potential and treat and simply 
deal with the continuous tail potential. The first term can be evaluated by introducing the 
inverse Laplace transform of rv ta u (r) , 

rvtau (r) = / ds exp (sr) V taU (s) (A4) 
Jo 

and likewise for rg hs (r) so that 

J d~f g hs (r) v taU (r) = 4tt^ dr r 2 g hs (r) v taU (r) (A5) 

roc roc 

= 4ir ds V ta ii (s) / dr rg hs (r) exp (-sr) 
Jo Jo 

= 4vr r dsVtau{s)G{s) 
Jo 

where G(s) is the Laplace transform of rg^ (r), which is known analytic function in the PY 
approximation 

G{s;d) = d 2 G PY {sd) (A6) 
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„ / \ xexp(— x) F (x 

Gpyix) 



Fix) 



1 + 12r/ exp (— x) F (x) 

1 1 + Ax 



127? 1 + (A - 1) x + (| - A) x 2 + (\A - ±gO x 3 



A 1 + V2 



l + 2r/ 

The integral in eq. ()A5|) is controlled by the exponential decay of G(s;d) and is easily per- 
formed numerically. Note that for the ten Wolde-Frenkel potential, we have that 



Vtail 0) 



960a 



2 



/ ( 'sof + 45 (so) 3 + (105 + 480a) so) sinh so 



V 



(A7) 



10 (so) 4 + (105 + 480a) (so) 2 ) cosh so 

The Percus-Yevick pair distribution function becomes exact at low densities but is only 
semi-quantitatively accurate at moderate to high densities. Compared to the pdf deter- 
mined from computer simulations, its oscillations are slightly out of phase and the pressure 
calculated from it is in error. The Verlet Weiss pair distribution function is a semi-empirical 
modification of the basic Percus-Yevick result designed to correct these flaws. It is written 
as 

9vw (r; p,d) = (r - d) (g PY (r; p, d ) + — exp (-m(r - d)) cos (m (r - d))) (A8) 

where the step function (r — d) ensures the fundamental property that the pdf vanishes 
inside the core, do is an effective hard-sphere diameter which has the effect of shifting the 
phase of the oscillations, and C and m are chosen to give the accurate Carnahan-Starling 
equation of state via both the pressure equation and the compressibility equation. To apply 
the Laplace technique in this case requires some care since what we know is the Laplace 
transform of gpy (r; p, d ) and not that of (r — d) gpy (r; p, d ). So we rewrite eq. (jA8|) as 

9vw (r; p, d) = g PY (r; p, d Q ) + (0 (r - d) - (r - d ))g PY (r; p, d ) (A9) 

C 

+0 (r — d) — exp (—m(r — d)) cos (m (r — d)) 
thus separating out the known PY contribution. This gives 

dPf g hs (r) vtaii (r) = J ' d~? g PY (r; p, d ) v ta u (r) (A10) 

rd 

+4tt / r 2 dr g PY (r; p, d ) v taiX (r) 

J do 

+47T / r 2 dr — exp (—m(r — d)) cos (m (r — d)) Vtau ( r ) 
Jd r 
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where the first integral can be evaluated via the Laplace transform technique, provided 
that do > a, the second integral is over a finite interval (for which one could analytically 
approximate the pair distribution function as in ref.j^J) while the third integral is easily 
evaluated numerically. All parts of the calculation are therefore well controlled. 

Finally, we note that the same techniques can be adapted to the evaluation of the second 
order contribution to the free energy. 
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